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1.  INTRODUCTION 


Whenever  a  beam  of  radiation  passes  through  an  aerosol  medium,  its 

energy  always  decreases  during  transit.  The  loss  of  energy  is  due  to 

either  scattering  or  absorption  or  both  in  the  medium.  The  term  extinction 

is  defined  as  the  sum  of  absorption  and  scattering.  The  transmittance  T 

for  a  parallel,  monochromatic  beam  of  radiation  of  wavelength  X  after 

passage  through  a  homogeneous  aerosol  layer  of  length  L  and  extinction 

coefficient  aext(^)'  in  the  single  scattering  (SS)  approximation,  is  given 

”T  ( A  ) 

by  Bouguer-Beer  law,  T  =  e  ,  where  the  optical  depth  t(A)  =  0  (A)  L, 

ext 

In  SS,  the  scattered  photon  reaching  the  detector  suffers  only  one 

scattering  event;  whereas,  in  multiple  scattering  (MS),  it  suffers  more  than 

th 

one  scattering  event.  Thus,  the  N  order  scattering  refers  to  a  case  when 

a  photon  is  scattered  and  rescattered  N  times  before  reaching  the  detector. 

The  basis  of  the  SS  theory  is  that  if  one  knows  the  scattering  properties 

of  an  individual  particle,  then  the  scattering  effect  of  N  similar 

particles  is  simply  N  times  that  of  a  single  particle.^"  Such  a  simple 

direct  proportionality  to  N  particles  is  what  makes  the  SS  theory  so  simple. 

In  the  case  of  MS,  this  proportionality  to  N  particles  no  longer  holds, 

thereby  making  the  theory  very  complicated.  Therefore,  except  for  a  few 
2-14 

researchers,  not  much  work  has  been  done  on  experimental  investigations 

on  MS;  even  though  relatively  more  work  has  been  done  on  the  theoretical 
15-75 

aspects  of  MS.  Therefore,  a  coordinated  research  program  was  undertaken 

between  the  Institute  for  Atmospheric  Optics  and  Remote  Sensing  (IFAORS)  and 
Army  Aberdeen  Chemical  System  Laboratory  (CSL)  to  implement  systematic 


7 


research  effort  in  the  theoretical  and  experimental  aspects  of  MS  in  aerosols 
in  controlled  laboratory  environment.  The  theoretical  research  was  performed 
at  the  Institute  and  the  experimental  work  at  CSL.  Important  features  of 
such  a  theoretical  and  experimental  approach  to  MS  investigations  were  that 
the  numerical  results  based  on  theoretical  MS  models  could  be  compared  with 
experimental  results  and  vice  versa.  Laboratory  experiments  can  serve  to 
model  real  atmospheric  conditions;  and  validation  of  MS  theories  are  essential 
prerequisites  for  gaining  confidence  in  their  applications  in  real  atmospheric 
conditions. 

In  this  program,  it  was  decided  to  start  with  the  simple  experiments  on 
measurement  of  attenuation  of  visible  and  IR  laser  beams  traversing  through 
fog  oil  aerosols  contained  in  a  closed  chamber,  such  as  making  measurement  of: 
(1)  optical  depth  for  three  or  four  wavelengths  as  a  function  of  time,  (2)  size 
distribution  of  aerosols,  (3)  intensity  of  scattered  radiance  as  a  function 
of  transverse  distance  from  the  direct  beam  axis. 

On  the  theoretical  aspects  of  the  program,  it  was  decided  to  develop 

22-24 , 47-48  ,  .  ^  63-69  ^  , 

exact  and  approximate  methods  based  on  solutions  of  the 

radiative  transfer  equation  in  the  small  angle  approximation.  Computer 

22-24  47-48 

programs  for  two  exact  methods  '  were  developed,  optimized,  and 

checked  for  consistency,  and  their  results  for  scattered  radiance  as  a 

function  of  optical  depth  for  the  various  orders  of  scattering  are  shown  in 

Section  2  of  this  report.  In  addition,  the  theory,  computer  codes,  and 

computations  for  radiance  as  a  function  of  the  detector  f ield-of-view  using 

63-69 

an  approximate  method  are  discussed  in  Section  3.  Conclusions  and 

recommendations  for  future  research  efforts  are  given  in  Section  4. 
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2.  EXACT  METHODS  FOR  COMPUTING  HIGHER  ORDER 
AND  TOTAL  MULTIPLE  SCATTERING  OF 

laser  beams 


2.1  INTRODUCTION 

Light  propagation  through  optically  thick  particulate  medium  is 
basically  a  multiple  scattering  problem  in  which  rays  or  photons  traverse  a 
medium  of  scatterers  and  undergo  many  scattering  events  before  escaping.  When 

a  collimated  laser  beam  is  incident  upon  a  thick  scattering  medium,  it 

16 

will  suffer  spatial  and  angular  broadening  as  well  as  attenuation.  These 

effects,  that  inhibit  optimum  performance  of  both  optical  and  communication 

2-14 

systems,  have  recently  become  a  focus  of  intensive  experimental  and 

15-75 

theoretical  studies.  References  17-24  devoted  to  analytical  investi- 

16  25  26 

gation,  parallel  the  Monte-Carlo  treatment  9  9  of  the  problem. 

The  natural  framework  to  study  the  multiple  scattering  problem  is 
the  transport  theory.  This  theory,  which  does  not  directly  include 
diffraction  effects,  deals  with  the  transport  of  energy  through  a  medium 
containing  particles.  Exact  general  solutions  to  the  equation  of  transfer 
have  not  been  obtained  to  date.  However,  there  are  some  special  cases 
where  simple  and  useful  approximate  solutions  are  available.  To  mention 
but  two  limiting  cases,  for  tenuous  particle  distribution,  the  first-order 
multiple  scattering  theory  can  be  used,  and  for  dense  distribution  the 
diffusion  approximation  is  appropriate.  When  particle  size  is  much 
greater  than  a  wavelength,  as  is  the  case  typically  for  fogs,  it  is  possible 
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to  simplify  the  equation  of  transfer  by  employing  the  small  angle  approxima¬ 
tion.  This,  in  particular,  involves  the  convolution  integral  in  the  gain 
term  (containing  the  scattering  phase  function),  which,  obviously,  lends 

itself  to  the  application  of  the  Fourier  analysis.  In  an  interesting 

2  8 

investigation,  Stotts  showed  that  the  use  of  Fourier  analysis  requires  the 
inclusion  of  all  position  and  angular  spatial  frequency  components  in  order 
to  ensure  the  correct  behavior  for  the  scattered  radiance  in  the  small-angle 
approximation.  His  results  preclude  the  low  order  truncation  in  the  series 
expansion  of  the  transformed  phase  function,  as  was  attempted  in  Refs.  18  and 
44.  The  high  frequency  limit,  which  is  related  to  the  first-order  scattering, 
would  not  be  adequate  as  the  optical  thickness  becomes  larger  and  larger. 
Therefore,  it  becomes  important  to  estimate  the  contributions  to  the  scattered 
radiance  coming  from  scattering  events  of  different  orders,  without  the 
limitations  imposed  by  the  approximate  Fourier  techniques. 

The  purpose  of  this  paper  is  to  provide  a  systematic  study  of  contri¬ 
butions  of  increasing  order  of  scattering  for  both  realistic  and  model 

aerosols.  Our  startina  point  will  be  the  theory  developed  by  Tam  and 
22-24 

Zardecki,  mainly  to  deal  with  multiple  scattering  of  laser  beams  in 

fogs.  This  theory  will  be  summarized  in  Section  2.2,  which  extends  the 
original  formulation  for  laser  beams  emanating  from  a  point  source  to  the 
more  realistic  case  of  laser  beams  of  finite  spatial  width.  Alternate  schemes 
for  solving  the  problem  are  described.  In  Section  2.3,  we  show  how  to 
construct  the  Gaussian  phase  functions  that  best  approximate  the  exact 
Mie  theory  phase  functions  for  both  monodisperse  and  polydisperse  aerosols. 

In  Section  2.4  we  discuss  the  computational  aspects;  and  in  Section  2.5  we 
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discuss  the  significance  of  the  computed  results  in  the  context 
of  laser  beam  propagation  in  dense  aerosols. 


2.2  SOLUTION  FOR  BEAM  IRRADIANCE 


Under  the  small-angle  approximation ,  and  with  the  z-axis  chosen  along  the 
direction  of  the  incident  beam,  the  equation  of  transfer  satisfied  by  the 
radiance  (or  specific  intensity)  distribution  function,  I  (^,£,_z)  ,  reads  as 


4) 


a  i 
e 


<J>')  I  ( r , z 1 ) 


dV 


(1) 


where  is  the  transverse  component  of  the  unit  propagation  vector;  the 
phase  function,  P,  is  considered  to  be  a  function  of  the  difference  of 
and  ' ;  and  0^  and  a  denote  the  extinction  and  scattering  coefficients, 
respectively.  In  terms  of  the  azimuthal  angle  p,  the  components  of  are 
(9  cos  p,  0  sin  p)  ,  where  9  is  the  scattering  angle.  In  addition,  we  have 
adopted  the  convention  in  which  the  location  of  a  given  point  in  space  is 
specified,  in  rectangular  coordinate  system,  by  (_r,z)  ,  where  £  =  xx  +  yy. 

Since  most  of  the  energy  scattered  by  large  particles  in  a  unit  volume 
is  scattered  in  the  forward  direction,  we  assume  that  in  the  forward  direction 
(or  forward  lobe)  the  angular  distribution  of  scattered  radiance  can  be  represented 
by  a  Gaussian  function.  Thus  the  phase  function  P(0)  may  be  written  in 
terms  of  the  projected  angles  <£  as  follows  (see  Fig.  1) : 

2  2 

P(£)  =  —  exp  (-a  $  )  (2) 


and  that  the  incident  laser  beam  can  adequately  be  represented  by  a 
Gaussian  functional  form,  both  for  the  spatial  distribution  and  the 
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z 


y 


sin 


FIG.  1.  The  scattering  geometrical  relationships.  In  the  small 
angle  approximately  the  components  of  <J>  are  (({)x,  <J>y)  = 
(8  cos  p,  0  sin  p) ,  9  being  the  scattering  angle ,  p 
the  azimuthal  angle,  n_  =  (({>_,  0)  ,  and  the  modulus  of 
equals  0,  i.e.,  |<j)J  =  0. 


6  sin<£,cos0) 
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angular  divergence.  Thus,  the  boundary  condition  at  z  =  0  is  expressed 
in  the  form 


=  l(<f>,r,0) 


p02  2  -2  ,  .2,2  2  2 

=  FB  y  it  exp  (-B  Q  -  Y  r  ) 


(3) 


where  F  is  the  total  power  received  at  z  =  0,  the  entrance  plane  of  the 
cloud  of  particles. 

Equation  (1)  may  be  solved  by  the  use  of  Fourier  transform  techniques. 
With  the  incident  beam  specified  by  Eq.  (3) ,  we  find  by  standard  methods20'29 
the  radiance  of  the  scattered  beam  in  the  form  of  a  Fourier  integral  over 

the  variables  (£,  n)  ,  which  are  the  Fourier  conjugates  of  variables  (r,^)  in 
real  space. 


I(<j>,£,z)  =  - T 

(2tt) 


S(E,n.z)  +  4*£  d2n 


(4) 


where 


I  (£,n , z)  =  F  exp(-Oez)exp 


(£  +  zn ) 


4  [" 


a 


exp  [f2(z)] 


(5) 


In  Eq.  (5) ,  the  function  fi(z)  is  related  to  the  Fourier  transform  of  the 
phase  function  through  the  relation: 


ft(z) 


z 

0  P  ( £  +  z  '  n)  dz ' 

5  Jo 


(6) 


which,  in  the  case  of  the  Gaussian  functional  form  given  by  Eq.  (2)  ,  becomes 


exp  < 

-  2  \L+  Hz‘\  } 

dz' 

(7) 

Jo  ' 

1  4a  J 
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As  has  been  shown  in  Ref.  22,  a  form  convenient  for  numerical  computation 
is  obtained  as  a  power  series  generated  by  expanding  the  exponential 
containing  fl(z)  in  a  Taylor  series  and  carrying  out  the  integration  over 
^  and  X)  before  z  1  . 

Because  of  the  peaked  nature  of  the  phase  function,  the  angular 

distribution  of  the  incident  beam  is  retained  through  several  scattering 
13 

events.  Therefore,  for  a  detector  whose  field  of  view  is  considerably 

larger  than  a  1,  which  represents  the  angular  spread  of  the  phase  function, 
the  irradiance  defined  as 


K$,r_,z)  d2<j> 


N(r,z) 


(8) 


will  effectively  be  equal  to  the  received  power.  In  the  following,  discussions 
we  shall  be  mainly  concerned  with  the  evaluation  of  the  beam  irradiance,  thereby 
reducing  the  number  of  independent  parameters  to  only  two,  namely,  r  and  z. 

On  combining  Eqs.  (4)  and  (8) ,  one  readily  obtains  the  irradiance 
in  the  form  of  the  following  two-dimensional  Fourier  integral 


N  ( r , z) 


(2tt) 


1 


(9) 


where 


I(ri,z)  =  I  (£  =  0,ri,z) 


(10) 


(11) 
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with 


z 

f  > 

n2  z'2 

aa  /  ir 

*»  -  »  —  S 

-  (D 

exp  ¥ 

.  2 

4a 

’  az  - 

n 

'0 

k  * 

(12) 


Eq.  (9)  can  be  simplified  by  the  following  two  approaches:  one,  gives 
the  separate  contributions  to  total  irradiance  due  to  both  single  and  higher  order 


scattering;  whence  the  second,  gives  the  total  multiple  scattering  (MS)  contri¬ 
bution.  It  should  be  noted  that  in  either  case  one  also  gets  the  unscattered 
contribution  to  total  irradiance. 


(i)  Various  Orders  of  Scattering  by  Series  Expansion 
On  expanding  the  exponential  term  containing  f2Q  ( 2 )  in  Eq.  (11)  and 
performing  the  integration  over  rj,  we  obtain 

N(r,z)  =  l  Nm(r,z)  *  (F/ttz2)  exp(~a  z)  l  *°s2'  N°(r,z)  (13) 

-  m=o  ~  e  m=Q — — m  - 

where  for  the  zeroth  order  (m  =  o)  scattering, 


o2  2  2 

N°(r,z)  =  — 2 


L..I tV  \, 

t  B2*vVr 


.2  2  2 


exp 


exp 


24-} 

W  (z)  > 


6  +  y  z  v  6“  +  y~z^'  w^z) 

With  W(z)//Tas  the  e"  2-point  Gaussian  width  at  axial  distance  z; 


for  the  nth  order  scattering 

1  1 
f  ( 

jO  , 


Nu(r,z)  = 
m  — 


dz. 


dz  A  (z)  exp 
mm 


(r/z) 


A  ( z 
m 


A 


(14) 


(15) 


with 


A  (z) 
m 


m 

I  (1-z.)2 

3-1  3 


and  for  the  single  scattering 


Nx (r , z )  = 


dz1A“1(z) 


+  6~2  +  Y~2z"2  (0  £  Zj  i  1) 

(m  =  1)  : 

exp [- (r/z) 2/A^  (z) ) 


(16) 


(17a) 
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with 
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\(z)  =  a-2  (1-z^) 2  +  6~2  +  y”2z"2 

It  should  be  noted  that  the  width  W 

o 

and  W(z)  by  the  relations 

Y  =  /2/W  ,  6  =  kW  //T~ 
o  o 

W2(z)  =  W2  +  (2z/kW  )2 
o  o 


(0  <  z,  i  1)  (17b) 

beam  waist  is  related  to  y,  & 

CL8a) 

(18b) 


The  problem  has  thus  been  reduced  to  a  numerical  computation  of  integrals 
given  by  Eq.  (15),  which  describes,  according  to  the  Wentzel  method,  ,23'30'71  the 
m-fold  scattering  contributions  to  the  observed  irradiance. 


Single  Order  Scattering  Contribution 

The  single  scattering  expression  obtained  from  Eqs.  (13,  15,  and  16)  by 

setting  m  =  1  is  difficult  to  evaluate  numerically  because  it  has  a  narrow  spike 


at  2j  1  Ec2'  0-6)  *  This  problem  can  be  solved  by  a  variable  transformation 


from  z  .  to  ,  namely 


2z4 


(1  -  z  . ) 
_ 3 _ 

2 

a 


w 


=  tan  6 


W2  (z) 


(19) 


which  for  single  scattering  (j  =  m  =  1)  leads  to: 

6 


N  (r , z)  =  N  1 (r , z) 
1  —  o  — 


~  i  e 

OJ _ 

AB 


— T 


exp 


2r  a2  •  2* 

— r-  A  sin  0 
W 


d6 


(20) 


where 


0) 


is  single  scattering  albedo; 


T  =  0eZ; 

6  -  tan  ^  AB  ? 

o 


and 


A  =  W  /W (z) 
o 

B  =  /2*  z/aw 

o 

N  (r ,z)  =  h°  (r,z) 


7TZ 


O  ~ 


16 


This  new  expression  in  Eq.  (20)  is  much  easier  to  evaluate  than  the  one  in 
Eq.  (17a,).  These  results  are  discussed  later. 


(ii)  Total  Multiple  Scattering  Contribution 

If  we  are  interested  in  obtaining  just  the  total  MS  contribution,  without 

the  need  to  know  the  individual  contributions  due  to  the  different  orders  of 

47  48 

scatterings,  then  we  could  use  an  alternate  formulation  '  obtained  by  simplifying 
Eq.  (9)  by  an  appeal  to  circular  symmetry  of  the  propagating  beam.  By  integrating 

the  two  dimensional  integral  of  Eq.  (9)  over  the  azimuthal  angle,  we  obtain  the 

-  n  .  .  25, 26 

following  expression: 


N(r,z) 


Kn,z)  jq  (nr)n 


dri 


(21) 


where  I  (r|,z)  is  defined  by  Eqs.  (10)  and  (11),  which  can  be  written  in  terms  of 
w  as 


o 


K  T\,z) 


(22) 


This  formulation  is  computationally  easier  to  evaluate  than  the  series  formula¬ 
tion  of  Eq.  (13) ,  provided  the  spatial  point  of  interest  lies  within  a  volume 
bounded  by  a  surface  coaxial  with  the  symmetry  axis  of  the  propagating  Gaussian 
beam  and  at  radial  distance  equal  to  the  width  of  the  unscattered  Gaussian  beam 
at  that  location. 


2.3  GAUSSIAN  PHASE  FUNCTION 

The  radiative  transfer  multiple  scattering  formulation  adopted  here  assumes 
that  the  scattering  phase  function  is  Gaussian.  In  general,  such  an  assumption 
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does  not  strictly  hold  true  for  scattering  profiles  obtained  by  applying  the 
Mie  theory  results  for  spherical  particles. 

Nevertheless,  an  assumption  of  Gaussian  phase  function  could  be  justified  in 
certain  domains  of  aerosol  particle  size  distributions.  When  the  Mie  scattering 
phase  function  is  highly  peaked  in  the  forward  direction,  then  the  intensity 
contributions,  resulting  from  scattering  a  few  degrees  from  the  forward,  are 
negligible.  Under  these  conditions  an  equivalent  Gaussian  phase  function  could 
be  fitted  to  the  forward  peak  of  the  Mie  scattering  phase  function. 

The  procedure  to  obtain  the  Gaussian  equivalent  phase  function  is  described 
as  follows. 


(i)  Exact  Mie  Phase  Function  Formulation 

It  is  assumed  here  that  the  individual  scatterers  are  isotropic,  homogeneous 
spheres,  which  enables  one  to  apply  the  results  of  the  Mie  theory  to  evaluate  the 
exact  form  of  the  volume  scattering  phase  function  P  (0  ,x,m)  [’sr~1]  ,  as  well  as 
ag(xfS)  [m  ]  and  Oe(x,m)  the  volume  scattering  and  extinction  coefficients, 
respectively,  where  0  is  the  scattering  angle  (Fig.  1)  ,  x  =  27Tr/A  is  the  size 
parameter,  r  £pm]  is  the  particle  radius,  A  [pm]  is  the  radiation  wavelength,  and 
m  =  m'  -  im"  is  the  aerosol  complex  refractive  index.  The  exact  scattering  phase 
function  is  defined  as 


P(0,x,m)  = 


27 TO 


f (0,x,m)  N  (for  monodisperse  aerosols) 

a 


(23a) 


where  0 

s 


2ttg 


s  r, 


=  Trr  Q  (x,rh)  N 
s  a 

f2  2 

=  Trr  Qs(x,m) 


2 

f (9,x,m)  n(a)da  (for  polydisperse  aerosols) 

(monodisperse) 
n(a)da  (polydisperse) 


and  f  (0,x,m)  =  (9,x,m)  +  (® m)^|  /2k 


(23b) 


(24a) 

(24b) 


(24c) 
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is  the  differential  scattering  cross-section 


f  2  -T1  . 

m  sr  I  ;  i 


and  i^  are  Mie 


intensity  functions;1  Qg  is  the  efficiency  factor;1  n(r)  is  the  aerosol  size 
distribution  | cm  ^  ym  1 


T  ~3 

cm 


C- 


;  and  is  the  monodisperse  aerosol  number  density 


?  and  r^  and  r^  are  limits  of  integration.  The  convention  introduced  in 
Eq.  (23b)  ensures  that  the  integral  of  P(0,x,m)  over  the  total  range  of  the 
argument  will  be  equal  to  unity. 

(iT)  Gaussian  Fit  to  Mie  Phase  Functions 


It  is  required  to  obtain  a  Gaussian  phase  function  of  the  form 

P(0,C t)  =  Ae" 


-e2a2 


(25) 


where  A  is  the  peak  value  of  the  phase  function  which  corresponds  to  the 

forward  scattering  peak  obtained  by  the  Mie  formulation;  and  a  is  the  inverse 

-1  2 
angular  width  (in  rad  )  of  the  equivalent  scattering  function,  a  is  obtained 

2  *  2 

by  plotting  the  normalized  Mie  phase  function  P(0)/P(O)  versus  0  deg  on  a  log- 
linear  scale  and  obtaining  the  slope  of  the  straight  line  that  best  fits  in  the 
least  squares  sense  the  forward  lobe  of  the  Mie  phase  function.  The  computations 
are  discussed  in  the  next  section. 

Such  an  approach  will  reproduce  the  optical  effects  of  the  forward  peak 
of  the  scattering  functions. 

2.4  COMPUTATIONAL  CONSIDERATIONS  AND  RESULTS 
(i)  Mie  Phase  Functions 

The  Mie  phase  functions  P(0,x,m)  were  computed  by  the  exact  Mie  scattering 

codes76  for  monodisperse  spherical  aerosols  of  various  radii  (in  the  range  2.0  ym 

77 

to  41.0  ym)  and  for  polydisperse  aerosols,  represented  by  analytic  models, 
such  as,  the  popular  Junge  power  law, 

n(a)  =  p^a  ^2  (a1  <_  a  (26) 
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76 

and  Deirmendjian  cloud  model  Cl;  here  P]_  and  p2  are  adjustable  constants  obtained 

from  observational  data.  These  normalized  Mie  phase  functions  are  plotted  as 
2 

functions  of  9  as  shown  in  Figs.  2-5  by  curves  or  straight  lines  connecting  the 
circles.  Figs.  2  and  3  are  for  the  case  of  monodisperse  aerosols?  and  Figs.  4 
and  5  are  for  polydisperse  aerosols.  In  the  case  of  particles  large  compared  to 
A,  the  Mie  phase  function  becomes  highly  oscillatory,  hence  the  dog-legged  shape 
of  the  plot  in  Fig.  3. 

(ii)  Gaussian  Fits  to  the  Forward  Lobe 

The  equivalent  Gaussian  phase  function  is  then  obtained  by  best  fitting  the 

functional  form  P(9,a)  in  Eq.  (25)  to  the  forward  lobe  of  the  normalized  Mie 

phase  functions.  The  best  fits,  obtained  by  the  use  of  a  least  squares  code, 

are  shown  by  the  dashed  line  in  Figs.  2-5.  The  slope  of  dashed  line  yields  a (deg  "S  ; 

and  the  normalization  constant,  the  value  for  A. 

Gaussian  phase  functions,  as  given  by  Eq.  (25),  were  best  fitted  to 

the  exact  Mie  phase  functions  for  monodisperse  aerosols  with  radii  in  the 

range  from  2.01  ym  to  40.2  ym  and  refractive  indices  m'  =  1.33  and  1.55. 

Essential  parameters  of  these  cases  are  summarized  in  Table  1.  The 

Gaussian  fit  obviously  supresses  the  oscillations  originating  from  the 

diffraction  effect.  These  oscillations  are  effectively  smoothed  out  in 

the  case  of  polydisperse  aerosols,  as  can  be  seen  from  Figs.  4  and  5. 

The  Gaussian  approximation  merely  leads  to  a  faster  convergence  for  large 

scattering  angles.  It  may  be  reminiscent  of  other  approaches  to  the 

45 

multiple  scattering,  e.g.,  the  statistical  theory  of  Whitney,  by  which 
the  fine  details  of  the  angular  scattering  problem  are  lost,  thus  circum¬ 
venting  the  need  for  the  generation  of  the  phase  function  with  a  Mie  code. 
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NORMALIZED  PHASE  FUNCTION  P(0)/P(O) 


FIG.  2 


Gaussian  fit  to  the  forward  lobe  of  the  normalized  Mie  phase  function 
for  X  =  0.425  mm  and  monodisperse  spherical  particles  of  a  =  10.0  pm 
and  m  =  (1.55  -  Oi) . 
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NORMALIZED  PHASE  FUNCTION  P(0)/P(O) 


FIG.  3.  Gaussian  fit  to  the  forward  lobe  of  the  normalized  Mie  phase  function 
for  A  =  0.45  ym  and  monodisperse  spherical  particles  of  a  =  40.2  ym 
and  m  =  (1.33  -  Oi) . 
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NORMALIZED  PHASE  FUNCTION  •'  P(0)/P(O) 


FIG.  4.  Gaussian  fit  to  the  forward  lobe  of  the  normalized  Mie  phase  function 
for  X  =  0.425  yim  and  polydisperse  aerosols  of  size  distribution  n(a) 
p^a“^  ,  with  p^  =  46.8,  0.4  ^  a  ^  13.0  ym,  and  m  =  (1.55  -  Oi)  . 
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FIG.  5.  Gaussian  fit  to  the  forward  lobe  of  the  normalized  Mie  phase  function 
values  obtained  from  (Table  35,  Ref.  21)  for  X  =0.45  pm,  Deirmend- 
j ian  Cloud  Model  Cl,  and  m  =  (1.34  =  Oi) . 
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TABLE  1.  Inverse  Angular  Width  a  of  Equivalent  Gaussian  Fit  to  the  Forward 


Lobe  of  the  Mie  Phase  Functions  for  Different  Monodisperse  Aerosols 


r  (|im) 

A  (um) 

X 

m' 

a 

deg"1 

rad*1 

2.01 

0.450 

28.06 

1.33 

0.29 

16.82 

5.00 

0.425 

73.92 

1.55 

0.71 

40.51 

7.75 

0.425 

114.58 

1.55 

1.00 

57.30 

10.00 

0.425 

147.84 

1.55 

1.35 

77.33 

15.50 

0.425 

229.15 

1.55 

2.15 

123.37 

20.00 

0.425 

295.68 

1.55 

2.60 

148.71 

26.80 

0.450 

374.20 

1.33 

3.58 

205.15 

31.00 

0.425 

458.30 

1.55 

4.12 

236.24 

38.00 

0.425 

561.30 

1.33 

5.05 

289.33 

40.20 

0.450 

561.30 

1.33 

5.05 

289.33 
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(iii)  Variation  of  a  as  Function  of  m'  and  x 

The  values  of  the  inverse  angular  widths  for  monodisperse  nonabsorbing 

particles  from  Table  1,  when  plotted  as  shown  in  Fig.  6,  fall,  as  expected, 

on  a  straight  line.  The  slope  of  the  lines  are  functions  of  the  index  of 

refraction  of  the  scattering  particles.  It  could  be  concluded  as  a  rule 

that  as  the  index  of  refraction  increases  for  nonabsorbing  particles,  the 

forward  scattering  peak  broadens  and  the  angular  width  per  degree  drops 

as  a  result.  Under  certain  conditions  it  is  difficult  to  estimate  accurately 

the  index  of  refraction  of  such  particles  or  it  may  be  an  unknown  quantity 

under  certain  circumstances.  On  examining  equations  (13)  through  (16)  we  can 

deduce  the  following  relationship 
~’N(a,z)  a  9n  (a ,  z) 

3cx  z  3z  (26) 


which  reveals  that  the  far  field  scattered  radiation  is  insensitive  to  small 
changes  in  cl,  the  inverse  angular  scattering  width. 

The  two  lines  in  Fig. 10  obey  the  following  fitted  expression: 


a 


3.25 


(m’) 


-1/2 


(27) 


where  a  is  the  inverse  angular  width  per  radian,  X  (wavelength)  and  r 

(particle  radius  in  pm;  and  m'  the  real  part  of  complex  aerosol  refractive  index). 

If  the  refractive  index  changes  by  a  small  amount,  then 


Act 


3.25 

-_r 


[•——  -  — — 1 


a 


Am1 

2m' 


which  means  that  the  percentage  change  of  the  inverse  angular  width  a  is 
half  that  of  the  corresponding  change  in  the  index  of  refraction. 


(28) 
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6 

w 

Cn 


(iv)  Computational  Problems  in  Evaluation  Eq.  (21). 

To  check  out  computation  scheme  for  solving  Eq.  (21),  we  decided  first 

to  compute  NQ(r,z)  by  setting  0^  =  0,  or  T  =  0,  in  Eq.  (22)  and  performing 

the  integrations  over  f|  for  various  values  of  r.  Results  for  N  thus  obtained 

o 

from  Eq.  (21)  using  a  1000-point  trapezoidal  integration  scheme  agreed  very 
well  with  those  from  Eq.  (14),  for  r  <  W(z),  i.e.,  the  discrepancies  were 
less  than  0.5%.  However,  negative  or  erroneous  results  were  obtained  for 
r  >  W(z) ,  due  to  the  need  to  include  a  very  large  number  (i.e.  over  20)  of 
cycles  in  the  evaluation  of  JQ(nr),  which  involves  values  of  Jq  comparable 
to  the  errors  in  the  numerical  representation  of  Jq  obtained  from  the  Handbook 
of  Mathematical  Functions,  p.  369  (Abramowitz  and  Stegun78) .  The  error  in 
the  use  of  this  representation  is  less  than  or  equal  to  1.6  x  10~8,  which 
implies  that  whenever  the  absolute  values  of  Jq  become  comparable  to  or  smaller 
than  this  error,  then  the  integration  of  Eq.  (21)  yields  erroneous  results. 

In  addition,  to  avoid  underflows  in  the  numerical  evaluation,  negative 
arguments  in  the  exponential  were  not  allowed  to  exceed  675.  As  a  result  of 
the  highly  complicated  numerical  integration  problems  arising  due  to  the 
oscillatory  nature  of  the  zero  order  Bessel  function  in  the  integrand  in  Eq.  (21)  , 
the  following  computational  method  was  developed  to  compute  the  total  irradiance 
N(r,z).  In  this  scheme,  first  we  determined  the  locations  of  the  zeros  of  the 
Bessel  functions,  and  the  intervals  between  the  zeros  were  then  subdivided  by 
a  predetermined  integral  value  to  accommodate  the  variations  in  the  Bessel 
function  cycles  within  a  fixed  range  of  integration.  This  method  is  very 
useful  due  to  the  fact  that  Eq.  (21)  becomes  highly  oscillatory  and  a  fixed 
interval  computational  scheme  gives  erroneous  results.  Computation  of  off- 
axis  irradiance  with  this  method  becomes  more  expensive  as  the  radial 
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distance  r  is  increased.  The  results  for  N(r,z)  are  given  in  Table  2. 


(v)  Computation  Problems  in  Evaluating  Eq.  (15)  for  Various  Scattering  Orders, 

(a)  First  Order  Scattering  Contribution 

Equations  (13) -(15)  were  found  to  have  a  spike  in  values  of  N  (r,z) 

m 

at  the  integration  limit  z^  =  1,  arising  due  to  the  expression 


C 


1 


1  +  (d  -  zx 


2  2 
Z  W 


(29) 


As  explained  in  Section  II (i)  ,  this  problem  was  solved  by  making  a 
variable  transformation  from  z^  to  5 ,  as  shown  in  Eq.  (19),  to  arrive  at 
transformed  Eq.  (20)  for  N^(r,z).  The  integral  in  Eq.  (20)  was  easily 
evaluated  by  using  a  simple  100-point  trapezoidal  integration  scheme.  The 
results  for  are  shown  in  Table  2. 


(b)  Second  and  Higher  Order  Scattering  Computations 

58  79 

For  higher  order  scattering  computations  the  Conroy  method  '  was 
adopted.  In  the  Conroy  method,  in  contrast  to  the  Monte  Carlo  method,  the 
sample  points  are  distributed  systematically  over  the  integration  region  and 
are  6ptimally  predetermined  to  achieve  convergence.  A  characteristic  feature 
of  the  Conroy  method  is  that  the  number  of  sample  points  will  need  to  be 
increased  to  achieve  a  specific  degree  of  accuracy  as  we  go  to  integrals  of 
higher  and  higher  multiplicity.  In  this  method,  eight  different  sample  point 
sets  are  available  for  multidimensional  integrations  for  m  =  2  to  12,  in  a 
tabular  form  in  Ref.  58  and  79. 

The  1154-point  sets  given  in  Ref.  58,  were  used  to  compute  2  to  7 
dimensional  integrations  (for  m  =  2  to  7) ;  the  3708-point  set  was  used  to 
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compute  8  dimensional  integration  (for  m  =  8) ;  and  the  3722-point  sets 
used  to  compute  from  9  to  12  dimensional  integrations  (for  m  =  9  to  12) .  The 
selection  of  these  point  sets  was  made  by  repeatedly  performing  the  integra¬ 
tions  with  the  different  sets  and  selecting  the  lowest  point  set  that 
achieved  convergence. 

2.5  DISCUSSION  OF  RESULTS  AND  CONCLUSIONS 

Table  2  presents  a  summary  of  the  results  obtained  from  Eqs.  (20),  (13), 

(15)  and  (21)  for  SS,  higher  order  scattering,  and  total  irradiance 

contributions  for  an  aerosol  cloud  represented  by  Gaussian  phase  function 
-1 

with  a  =  46.8  rad  ,  shown  in  Fig.  5,  at  transverse  distances  r  =  0.0  and 
2.0  cm  and  layer  depths  z  =  10~4,  10~3,  .  .  .,  102  and  103  km. 

In  Table  3,  we  show  a  typical  set  of  results  obtained  from  Eqs.  14  and  15 
for  the  irradiance  contributions  due  to  the  different  orders  of  scattering  for 
a  10  m  deep  cloud  of  aerosols  represented  by  a  Gaussian  phase  function  with 
T  =  289.33  rad  ,  shown  in  Fig.  3,  and  for  optical  depths  T  ■  0.1,  0.3,  1.0,  3.0, 
5.0  and  7.0.  Also,  shown  are  results  for  the  total  scattered  and  the  direct 
attenuated  irradiances.  Figures  7  and  8  depict  the  unscattered  (NQ) ,  MS  and 

total  scattered  irradiance  as  a  function  of  the  relative  transverse  distance, 
r/WQ  at  z  =  10  m  for  an  aerosol  cloud,  represented  by  a  Gaussian  phase  function 
with  a  =  289.33  rad  1  in  Fig.  3,  and  having  optical  depths  T  =  0.1  and  2.0, 
respectively.  Figure  9  depicts  the  same  quantities  for  aerosol  cloud  with 
a  =  33.07  rad  and  T  =  0.2  and  4.2.  Figure  10  shows  the  variations  of 
unscattered,  multiple  scattering  orders  (m  =  1,  2,  .  .  .  15),  total  MS,  and 
total  irradiance  received  by  an  open  detector  placed  coaxially  along  the  direction 
of  beam  propagation  for  various  optical  depths  up  to  j  =  10.0  due  to  an  aerosol 
cloud  having  a  =  289.33  rad  3  and  z  =  10  m. 
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Scattered  Irradiance  at  Various  z  and  r  for  Gaussian  Laser  Beam  of  Width  W  =  1.0  cm  Propagating  Through 
Water  Cloud  Cl,  Represented  by  a  Gaussian  Phase  Function  with  a  =  °46.8  rad--*-  shown  in 
Fig.  5  and  Having  an  Optical  Depth  T  =  1.0 
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Arguments  in  parenthesis  denote  powers  of  10. 
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RELATIVE  IRRA  Dl  ANCE  f  W  cm 


FI', . 


1  /  v- J/r^^^ince  as  a  function  of  the  normalized  transverse  dis 
,rAlo>  '"r  =  .10  289.33  rad'1  and  T  •  0.1. 
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RELATIVE  IRRADIANCE ,  W  cm 


FIG.  8.  Relative  irradiance  as  a  function  of  the  normalized  transverse  dis¬ 
tance  (r/WQ)  for  z  =  10  m,  a  =  289.33  rad-1  and  T  =  2.0. 
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NORMALIZED  IRRADIANCE 


FIG. 


.  Normalized  irradiance  as  a  function  of  the  normalized  transverse 
distance  (r/WQ)  for  z  =  10  m,  a  =  33.07  rad-1  for  T  =  0.2  and  4.2. 
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IRRAOI ANCE  ,  Wcm  pm 


Irradiance  vs.  T 


z  =  10  m 
tt  =  289  rad 


FIG.  10.  The  relative  irradiance  as  a  function  of  the  optical  depth  T, 

showing  the  unscattered  contribution  (denoted  by  0) ,  contribution 
due  to  various  orders  of  scattering  (denoted  by  1,  2,  3  ...),  total 
scattered  contribution  (blank  circles)  and  total  scattered  plus 
unscattered  contributions  (solid  circles) . 
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In  this  section  we  show  the  importance  of  higher  order  scatterings  as 
related  to  the  transmission  of  Gaussian  laser  beams  through  dense  scattering 
media.  The  Conroy  scheme  used  to  obtain  results  from  Eq.  14  differs  from  the 
Monte  Carlo  approach  in  that  the  points  of  integration  are  predetermined  to 
achieve  best  convergence  so  that  convergence  is  achieved  at  a  faster  rate  than 
the  Monte  Carlo  method,  resulting  in  a  lesser  computational  time.  The  formula¬ 
tions  described  earlier  are  based  on  the  Gaussian  representation  of  the 
scattering  phase  function  in  the  forward  direction,  in  which  the  essential  fac¬ 
tor  is  the  width  of  the  forward  scattering  lobe  for  monodisperse  or  polydisperse 
scatterers.  Using  the  Conroy  numerical  scheme  we  have  been  able  to  obtain  good  * 
theoretical  results  for  optical  depths!  £7.0.  At  T  =7.0,  the  contribution 
from  the  sixth  order  scattering  peaks  while  at  the  same  time  scattering  orders 
of  up  to  12  contribute  appreciably  to  the  total  scattering  (Fig.  10) .  Compu¬ 
tations  by  the  series  method  for  the  off-axis  locations  (e.g.,  r  =  2.0  cm)  of 
the  open  detector  require  the  same  number  of  integration  points  to  achieve  the 
same  accuracy  as  required  for  the  on-axis  (r/W^  =  0)  case;  whereas,  by  the 
Bessel  function  integration  method,  the  number  of  integration  points  need 
to  by  increased  for  the  off-axis  locations,  as  explained  earlier. 

On  the  other  hand,  the  alternative  formulation  of  Eq.  (21),  involving 
Bessel  functions,  cannot  estimate  the  contributions  of  the  different  orders  of 
scatterings  although  it  is  a  powerful  tool  because  it  can  compute  the  total 
scattering  irrespective  of  the  value  of  the  optical  depth  as  long  as  we  are  within 
,a  volume  bounded  by  the  surface  described  in  Section  2.2  (ii) . 
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From  Fig.  10,  we  notice  that  for  x  0.1,  only  SS  contributes  to  the 
total  irradiance;  for  0.1  X  _<  0.3,  the  second  order  contributions  also  needs 
to  be  included;  and  for  T  >  0.3,  higher  orders  of  scattering  need  to  be  taken 
into  account.  For  example,  for  T  -  2.0,  the  maximum  contribution  to  irradiance  is 
made  by  the  first  order  scattering  and  that  contributions  from  up  to  5  orders  of 
scattering  are  significant.  Similarly,  for  T  ~  7.0,  the  maximum  contribution  is  made 
by  the  6th  order  scattering,  and  that  contributions  from  up  to  13  orders  of  scatter¬ 
ing  are  significant.  The  contributions  from  the  various  orders  of  scattering  seem 
to  be  enveloped  by  a  curve  which  has  a  mode  at  about  T  -  1.0,  and  falls  monotomically 
for  T  >  1.0.  Similarly,  the  total  scattered  irradiance  is  given  by  a  curve  similar 
of  similar  shape,  which  as  a  mode  at  about  T  —  1.0  and  falls  monotomically  for 
T  >  1.0  but  at  a  lesser  rate  than  the  previous  envelope.  The  total  irradiance 
received  by  the  detector  is  given  by  the  curve,  obtained  by  adding  the  unscattered 
and  scattered  irradiances.  The  detector  receives  forward  scattered  radiation 
in  addition  to  the  direct  radiation,  thereby  causing  deviations  from  the  Beer- 
Bouger's  law  (i.e.,  transmittance  T  *  exp  (-i)),  which  is  represented  by  the 
dash  line  marked  0.  This  shows  succinctly  that  the  optical  depth  obtained  from 
the  transmissometer/radiometer  measurements  would  not  be  the  true  value  but  an 
apparent  optical  depth.  Therefore,  in  order  to  obtain  the  true  optical  depths , 
a  forward  scattering  correction  must  be  made  to  the  transmissometer/radiometer 
measurements.  A  detailed  discussion  of  the  MS  corrections  as  functions  of -the 
detector  field  of  view  and  optical  depths  is  given  in  Section  3  next. 
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3.  MULTIPLE  SCATTERING  CORRECTIONS  AS 


FUNCTIONS  OF  DETECTOR  FIELD  OF  VIEW 


3 . 1  INTRODUCTION 


In  making  optical  extinction  measurements  by  laser  beams  traversing 

aerosols,  an  error  is  always  made  because  some  forward  scattered  radiation 

invariably  enters  the  detector's  finite  field  of  view.  To  determine  the 

correction  due  to  multiple  scattering  requires  the  solution  of  the  radiation 

transfer  equation.  A  fruitful  approach  to  the  problem  has  been  a  Monte  Carlo 

technique,  because  it  can  handle  strange  geometries  as  well  as  inhomogeneities^ ' 

Recently,  considerable  attention  has  been  paid  to  analytical  solutions  to  the 

equation  of  transfer  in  a  form  appropriate  for  the  laser  beam  propagation 
..  22,24,47,48,60,61  ,  , 

problem.  Although  exact  solutions  have  not  been  obtained  to 


date,  there  are  some  special  cases  where  simple  and  useful  approximate  solutions 
to  the  equation  of  transfer  are  available.  For  tenuous  distribution  of  scatterers, 
the  first-order  multiple  scattering  theory  can  be  used,  and  for  dense  distribution 
the  diffusion  approximation  is  appropriate.  If  the  particle  size  is  large 
compared  with  incident  wavelength,  the  energy  scattered  by  the  particle  is 
largely  confined  within  a  small  angle  in  the  forward  direction  and,  therefore, 
by  employing  the  small-angle  approximation  it  is  possible  to  simplify  the 
equation  of  transfer.  In  Ref.  61,  a  systematic  study  of  contributions  of 
increasing  order  of  scattering  for  both  realistic  and  model  aerosols  has  been 
conducted. 
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27  62 

As  indicated  by  Ishimaru,  '  the  first-order  multiple  scattering 
approximation  is  applicable  when  the  density  of  scatterers  is  so  low  that 
the  diffuse  (incoherent)  intensity  is  considerably  smaller  than  the  reduced 
(coherent)  intensity.  This  will  certainly  be  the  case  if  optical  distance 
traversed  by  the  bean  is  much  smaller  than  unity.  However,  the  same  weak 
fluctuation  case  is  also  encountered  in  the  situations  where  the  receiver 
has  a  narrow  receiving  angle.  In  this  case,  the  amount  of  scattered 
intensity  entering  into  the  receiver  is  small  compared  with  the  direct 
coherent  intensity,  and  therefore  the  received  field  is  predominantly 
coherent.  The  effect  of  the  detector's  finite  field  of  view  on  the 
received  power  has  been  a  subject  of  comprehensive  investigations  related 

O  p  p  OA 

to  forward  scattering  corrections  for  optical  extinction  measurements. 

The  purpose  of  this  effort  was  to  study  the  effect  of  a  finite  field  of 

view  on  the  intensity  and  the  received  power  of  a  laser  beam  undergoing 

multiple  scattering.  Our  analytic  approach  is  based  on  the  theory  of 

Dolxn  and  Fante,  summarized  m  Section  3.2.  In  Sections  3.3  and  3.4, 

we  apply  this  theory  to  Gaussian  beams.  Numerical  results  relevant  to 
the  beam  propagation  in  a  water  cloud  and  model  aerosol  particles  are 
presented  in  Section  3'.  5. 

3.2  THE  FANTE-DOLIN  THEORY 


Our  considerations  will  be  based  on  the, equation  of  radiative  trans¬ 
fer  for  the  radiance  (specific  intensity)  distribution  function,  I((£,r,z). 
In  the  small  angle  approximation,  I(cj>,r,z)  satisfies 


1  9l  .  3l 

<j>  *  +  -r—  +  01  =  0 

or  dz  s 


F($  -  I' )  !($' ,r,z)  d 


(la) 
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where  o  and  G are  volume  extinction  and  scattering  coefficients  (m  for 
the  aerosol  medium;  r  is  the  component  of  the  position  vector  transverse  to 
the  z  axis;  and  is  the  transverse  component  of  the  unit  propagation  vector 
as  shown  in  Fig.  1  in  Section  2.2. 

In  the  work  of  Dolin64  and  Fante, 20,65,66,69  which  applies  to  sharply  peaked 

-V  •> 

phase  functions,  I(<J>',r,z)  in  the  integrand  on  the  right  hand  side  of  Eq. 
is  expanded  in  Taylor  series  about  <£'  =  <£.  More  precisely,  I(<{i,r,z)  is 
split  into 


T  frh 


(30) 


where  superscripts  o  and  s  refer  to  unscattered  and  scattered  radiance, 
respectively.  Then,  one  expands 

I(S*  (<}>,r,z)  =  ($,r,z) 

+  Z  (<j>*  -  <{>)  1^  (<£,r,z) 

k  k  d\ 

+  1  llg  ~  ^k  (<J>'  "  l(S)  (^,z)  +  *  *  *  (31) 

x  36  k  Jc 


where  k  and  £  refer  to  the  Cartesian  components  of  the  <j>  vector.  Recall 


that  <j>  =  <f>  x  +  <j)  y. 

x  y 


Equation  (1)  now  yields 


.(o) 


3 1 


(o) 


9  z 


+  a  i 


(o) 


=  0 


(32) 
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while  I 


( s )  . 

is  determined  from  the  nonhomogeneous  equation 

-  ?  <♦’»  vj  i“>  ♦ 


.  (s) 


+  a(l  -  U)  l'°'  =  00)  I  p(<jt)  -  $•)  l(o)  ($'  ,r,z)  d^<j>'  (33) 


where  w  =  a  /a  is  the  single  scattering  albedo;  and  <|2>  is  defined  as 


42>  = 


P(l)  %2  d2  4> 


(34) 


In  deriving  Eq.(33)  we  have  neglected  the  terms  of  higher  order  than  second 
in  Eq.  (31)..  The  solutions  to  Eqs.(32)  and  (33)  can  be  obtained  by  the  method 
of  characteristics  in  the  Fourier  space  of  the  variables  $  and  r.  With 
the  Fourier  transforms  defined  generically  as 

CO 

a2*  d2r  I(|,;,z)  e11^  * 


K5,n,z)  = 


(34) 


P(Q  = 


d2<J>  P(|)  eU** 


(35) 


we  obtain  the  following  solutions 


*(o)  it  t  \ 
i  (t>.-n,z) 


=  i(5  +  nz,n,z  =  o) 


-OZ 


(36) 


and 


*  (s)  -t  -*■  » 

i  (C»n»z)  =  i(5  +  nz,n,z  ~  o) 


P[|  +  n(z-z’ ) ]  e 


-oz 


•  exp 


z  1 


r00)  -f2  i±  1 2 

[“7“  <(i>  >  |C  +  n  (z  -  zn)  |  +  o(i-u))]  dzu 


dz  1 


(37) 


for  the  unscattered  and  scattered  contributions,  respectively. 
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3.3  CASE  OF  COLLIMATED  GAUSSIAN  BEAM 

We  assume  now  that  in  the  single  scattering  theory  the  scattering  phase 
~y 

function  P(4>)  is  given  by  a  Gaussian  function,  i.e.. 


P(^)  =  ~~  exp  (-a2  4 )2) 


(2) 


and  that  the  incident  collimated  beam,  directed  along  the  z  axis,  has  a 

Gaussian  spatial  form,  i.e., 

T/t  ■+  ..  _  -1  2  x  (2)  ,1,  ,  2  +2. 

I(q>,r,z  =  0)  =  FqTT  y  6  (<p)  exp  (-y  r  )  (38) 

(2)  . 

where  6.  is  the  Dirac  delta  function. 

In  the  situation  modeled  by  Eqs.(2)  and  (38),  we  obtain,  after  performing 

the  inverse  Fourier  transforms  of  Eqs.  (36)  and  (37) ,  the  explicit  formulas 

(o)  (s)  “►  *4- 

for  I  (<J>,r,z)  and  I  (<p,r,z).  They  read 


T  ^°)  I~X  Z  \  p  ~az  -Y2  X2  jp  (2) 

I  (<p,r,z)  =  F  ir  e  e  '  6 

o 


(I) 


and 


,  .  ^  __  F  awe 

I  (<p,r,z)  =  - 

(2-n) 


-az 


-1 


dz'  [4  A(z’)  C(z')  -B  (z ' ) ] 


(39) 


•  exp 


-*2  *4-  ->-2 

A(z')r  -  Bjz'l^r  +  0(2*  )0 

4A ( z ' ) C ( z ' )  -  B2(z') 


ujoz  ' 


(40) 


where  the  functions  A,  B,  and  C  are  defined  as 

1  +  awz’ 


A(z')  = 


B(z')  = 


4a 


2z’  +  awz' ' 


4a 


C(z' )  ,  -1^  + 

4y  4a 


(41) 
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(42) 


In  polar  coordinates  (0,b)  ,  the  vector  0  can  be  expressed  by 


0  =  9  cos  b 

x 


0=0  sin  b 

y 

For  a  detector  having  the  field  of  view  (FOV)  half-angle,  0  ,  the 

received  intensity  is  obtained  from  Eq.  (.40)  after  integration  over  the 

angle  b  in  the  range  (0,2tt)  and  over  0  in  the  range  (O#0  ).  The  received 
( s )  -> 

intensity  F  (0p,r,0z)  corresponding  to  the  scattered  beam  thus  becomes 


er 


( s)  -*• 

F  .  (0D,r,az)  = 


I(S) (0,r,az) e  d0 


(43) 


whe  re 


/  \  F  0)  e 

T  ( s)  / q  N  o 

I  (0,r,az)  =  - — 

2tt 


-a  z 


rOz 


[4A(|)C(f) 


B2(f>] 


-l 


Xo  [' 


0r  B (x/g) 


4A(x/a)C(x/a)  -  B  (x/0) 


r  A(x/a)r2  +  c(x/a)  e2  , 
exp  [ - - ] 

4A(x/a)c(x/a)  -  b  (x/a) 


cox  , 

e  dx 


(44) 


Here  Io  is  the  Bessel  function  of  zeroth  order,  and  x  is  a  dimensionless 

integration  variable. 


The  contribution  coming  from  the  unscattered  part  of  the  beam  is 

2+2 

(o)  /ft  -*■  -1  -oz  -Y  r 

F  (d^irtoz)  =  F  7T  e  e 
D  o 


(.45) 


It  is  thus  seen  that  the  numerical  results  can  be  extracted  from  this 
theory  in  a  rather  simple  manner  by  performing  merely  a  double  integration. 
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3.4  CASE  OF  GENERAL  GAUSSIAN  BEAM 

More  generally,  one  can  account  both  for.,  the  spatial  and  angular  divergence 
of  the  laser  beam  assuming  a  Gaussian  law  of  the  form 


,2  2 


I(<t>,r,z  -  0)  =  Fq  8  Y*  Tr  z  exp (-  <t>2  -  y2  r2)  (46j 

for  the  incident  beam.  Here,  the  parameter  B  describes  the  angular  divergence. 

By  following  exactly  the  same  line  of  reasoning  as  in  Section  III,  we 
obtain  the  following  contributions  F (s) (0D,?,az)  corresponding  to  the  scat- 
tered  beam,  viz., 


F(s> (0D,Yr,az)  =  0D 


I(S)  (0'0D,yr,az)  0'  d0 ' 


(47) 


where 


I(s)  (0,yr,az)  =  F  uoze  °z/2tt 

O 


x  |  [4A((x/a>  za}  c { (x/a)  za}  -  b 2{  (x/a)  za}j 


-l 


X  i 


0(yr)  -  B ((x/a)  za)/Y 


4a{ (x/a)  za}  c{ (x/a)  za}  -  b2{ (x/a)  za} 


x  exp 


(yr)2  a ((x/a)  za)/y2  +  c((x/a)  zo)62 

4a{ (x/a)  za}  c{ (x/a)  za}  -  b2{ (x/a)  za} 


x  exp(a)  xa  z)  •  dx' 

and  IQ  is  the  modified  Bessel  function  of  zeroth  order. 


(48) 
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For  the  unscattered  bean,  we  obtain 


,(o) 


(9D,Yr,ozj  = 


2F  8  exp(-Oz) 


eD2A 


rl 


6’  <30 '  lQ{2(y/a)  0^0 '  (Yr)  Oz} 


+  exp  {-  (Yr) 


(Y/a)2  (az)2  +  82 


eD2  e-2} 


(49) 


The  integrals  in  Eqs.  (48)  and  (49)  have  limits  0  and  1.  This  choice  of 


limits  anticipates  the  use  of  Gaussian  quadrature  in  numerical  computation. 


3.5  NUMERICAL  RESULTS 

The  numerical  values  for  the  on-axis  intensities  are  given  in  Tables  4-7 
and  the  corresponding  plots  are  given  in  Figs.  11-14,  in  which  we  show  the 
unscattered  (reduced)  intensity  marked  by  squares,  scattered  (diffuse) 
intensity  marked  by  circles,  and  total  intensity,  by  triangles,  on  the  beam 
axis  for  optical  depth  T  =  1,  4,  7,  and  10,  respectively.  A  divergent 
beam  with  the  parameters  6  =  277/ (YA)  and  Y  =  1.0  cm-1  is  assumed  to  be 
propagating  in  the  Deirmendjian  water  cloud  (Ref.  27)  model  Cl.  For 
A  =  0.45  ym,  the  extinction  coefficient  is  obtained  from  the  Mie  theory 
computations  when  a  modified  gamma  distribution  is  taken  for  water  cloud 
particle  size  distribution.  One  also  obtained  a  =  46.80  rad-1  for  the 
Gaussian  fit  to  the  phase  function  as  shown  in  Fig.  5.  For  a  coaxial 
detector,  with  a  diameter  RQ  =  1.0  cm,  the  same  as  the  laser  beam  diameter, 
the  intensity  F  becomes  a  function  of  Qp  and  z.  We  normalize  the 
intensities  by  dividing  out  the  factor  F^  (®D=00»  z  =  0)  .  For  the  sake 
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ON- AX  IS  INTENSITY 


Fig.  11:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 

corresponding  to  Table  4  for  a  =  46.80  rad-1,  RD  =  1.0  cm,  and  T  =  1.0. 
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ON-AXIS  INTENSITY 


Fig.  12:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 

corresponding  to  Table  5  for  a  =  46.80  rad"1,  rd  =  1.0  cm,  and  T  =  4.0 
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ON- AX IS  INTENSITY 


Fig.  13:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 

corresponding  to  Table  6  for  a  =  46.80  rad"1,  RD  =  l.o  cm,  and  T  =  7.0. 
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ON-AXIS  INTENSITY 


^^-9-  14:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 

corresponding  to  Table  7  for  a  =  46.80  rad-1,  RD  =  1.0  cm,  and  x  =  10.0. 
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of  clarity  in  graphical  presentation,  we  also  scale  the  normalized 
intensities  by  the  multiplicative  factor  exp(l). 

In  order  to  investigate  the  significance  of  the  optical  thickness,  we 
calculated  the  power  received  by  a  coaxial  detector  of  radius  =  1  cm. 
if  F(9D'r'T)  denotes  generically  the  beam  intensity,  then  the  received 
power  is  defined  as 

P  (Rd)  =  2TT 

The  received  power  is  also  scaled  by  the  multiplicative  factor  exp  (T) , 
and  the  power  received  at  z  =  0  is  divided  out.  The  numerical  values  for 
received  power  for  a  coaxial  detector  of  diameter  RD  =  1  cm  are  given  in 
Tables  4-7,  and  the  situation  is  depicted  in  Figs.  15-18  corresponding  to  the 
same  parameters  as  in  Figs.  11-14,  respectively.  As  the  FOV  increases,  both 
the  intensity  and  the  received  power  saturate  rapidly,  independently  of  the 
optical  depth.  This  corresponds  to  the  situation  of  an  open  detector.  The 
contribution  of  the  scattered  power  becomes  dominant  in  the  saturation 
region  for  optical  depth  of  the  order  of  10. 

Gaussian  phase  functions  as  given  by  Eq.  (2)  were  best  fitted  to  the 
exact  Mie  phase  functions  for  monodisperse  aerosols  with  radii  in  the  range 
from  2.01  to  40.2  pm.  In  Tables  8-11  and  Figs.  19-22,  we  show  the  intensity 
vs.  detector's  FOV,  and  in  Tables  8-11  and  Figs.  23-26,  the  received  power 
vs.  FOV.  Both  these  sets  of  results  were  computed  for  the  model  particles 

-1  _2  -i 

characterized  by  the  parameters  a  =  289.33  rad  ,  and  O  =  5.0  x  10  cm 


F(0D,r,T)  r  dr 


(50) 
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FOV  (RAD) 


Fig.  15:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of  detector 
FOV  corresponding  to  Table  4  for  a  =  46.80  rad~l,  RD  =  1.0  cm,  and 
T  =1.0. 
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RECEIVED  POWER 


Fig.  .16:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of  detector 
FOV  corresponding  to  Table  5  for  a  =  46.80  rad~l,  RD  =  1  0  cm  and 
T  =  4.0. 
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RECEIVED  POWER 


Fig.  18:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of  detector 
FOV  corresponding  to  Table  7  for  a  =  46.80  rad-1,  Rn  =  1.0  cm  and 
T  =  10.0. 
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TABLE  8.  Received  Power  and  On-Axis  Intensity  as  Functions  of  the  Detector  FOV  for  Large  Size 
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optical  depth=  10.00 
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10+ 


ON- AX IS  INTENSITY 


Fig.  19:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 

corresponding  to  Table  8  for  a  =  289.33  rad-1,  RD  =  1.0  cm,  and  T  =  1.0. 
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ON- AX  IS  INTENSITY 
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Fig.  20:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 

corresponding  to  Table  9  for  a  =  289.33  rad-1,  RD  =  1.0  cm,  and  T  =  4.0. 
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ON- AX  IS  INTENSITY 


Fig.  21:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 
corresponding  to  Table  10  for  a  =  289.33  rad“l,  RD  =  1.0  cm,  and 
T  =  7.0. 
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ON-AXIS  INTENSITY 


Fig.  22:  Normalized  intensity  on  the  beam  axis  as  function  of  detector  FOV 
corresponding  to  Table  11  for  a  =  289.33  rad-1,  RD  =  1.0  cm,  and  . 
T  =  10.0. 
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Fig.  23:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of 

detector  FOV  corresponding  to  Table  8  for  a  =  289.33  rad-1,  RD  =  l.o 
cm,  and  T  =  1.0.  Received  power  is  multiplied  by  exp  (T) . 
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RECEIVED  POWER 


Fig,  24:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of 

detector  FOV  corresponding  to  Table  9  for  a  =  289.33  rad-1,  Rq  =  1.0 
cm,  and  x  =  4.0.  Received  power  is  multiplied  by  exp  (x) . 
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RECEIVED  POWER 


Fig.  25:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of 
detector  FOV  corresponding  to  Table  10  for  a  =  289.33  rad~l,  Rq  = 
1.0  cm,  and  T  =  7.0.  Received  power  is  multiplied  by  exp  (t) . 
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RECEIVED  POWER 
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Fig.  26:  Received  power  (arbitrary  units)  on  the  beam  axis  as  function  of 
detector  FOV  corresponding  to  Table  11  for  a  =  289.33  rad-1,  Rjr,  = 
1.0  cm,  and  T  =  10.0.  Received  power  is  multiplied  by  exp  (t) . 


70 


at  X  =  0.45  ym,  for  a  coaxial  detector  of  diameter  R^  =  1  cm.  Note  that 
the  phase  function  and  the  Gaussian  fit  to  the  forward  lobe  of  the  Mie 
phase  function  are  plotted  in  Fig.  3. 

In  order  to  study  the  effect  of  the  detector’s  radius  R^,  we  let 
R D  assume  values  smaller  and  larger  than  the  beam  diameter  of  the  1.0  cm, 
viz,  R^  =  0.2  cm  and  2.0  cm.  Results  were  obtained  for  same  parameters  as 
in  Tables  8-11.  However,  for  the  sake  of  clarity  we  present  only  the 
results  for  T  =  4.0  for  the  two  cases,  viz,  R^  =  0.2  cm  (Table  12, 

Figs.  27-28)and  R^  =  2.0  cm  (Table  13,  Figs.  29-30).  Table  12  and 
Figs.  27-28  depict  the  corresponding  values  of  intensity  and  received  power 
as  functions  of  FOV  for  RQ  =  0.2  cm,  whereas.  Table  13  and  Figs.  29-30 
depict  the  corresponding  values  for  R^  =  2.0  cm.  The  three  cases 
(Rp  =  1.0,  0.2  and  2.0  cm;  T  =  4.0)  for  large  particles  with  phase  function 
parameter  a  =  289.33  cm  ^  show  clearly  the  contrast  in  the  behavior  of  the 
intensity  and  the  received  power. 

The  Fante-Dolin  approximation  employed  in  this  paper  enables  one  to 
estimate  the  corrections  to  the  Bouguer-Beer  law  for  a  receiver  having  a 
finite  field  of  view,  as  for  example  in  Ref.  45. 
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TABLE  12.  Received  Power  and  On-Axis  Intensity  as  Functions  of  the  Detector  FOV  for  the 


Onh-<Nr^r)°ocN<x)cy>'»-n^j-rfininininu)c£)£SiDiotp 

^’-^cMCNnn^Tr'cinininininininintninininiriinin 


o 

II 

h* 


S008888888888888888888888 

♦  ▲  ▲  x  x  ▲  4  ▲ 

£rj2JJi£?™S£!C£ECN00r')T"oo^ino,,"£CN0)n 

^cNrLnoi^rro3nninOT*conu5coaiOOO^»-»-’- 

oinincMc-cooorMcnai-'-n^^rinininininintniniptDip 

C7)CY)r-^-*-cNcNnnnTr«g-TrTfT3-'tf*ij-^^-^tTf53-xr'^-^ 


Tl 

C 

aJ 


CN 

o 

ll 

aP 


'  1  1  *  1  1  •  1  1  I  1  1  )  p  1  l  I  I  1  (  i  I  1  I  1 

/£,2i!J3!2jS!®®®®©©©©©®©©©©©ffi©©© 

^®®®®®tO®tP®U>tD(0C£)tDtO<X)(£)(£)tOlOtf>t£)y3t£J 

2?2£2?2?5?®®®aooo®®®®®®®®®®®®oooo 

®C7>Ol®C>CT>O>O>cn(T>O^O>®CT>0>C>Cy»OTO^O^®C7»O>OTO^ 

CP®(r>cj)oi®crcT»ocr>o,>o^a>oia)o^cT)a»cj)CT>o'a)oio^® 


U 

o 

-p 

ft 

a> 

o 

X 

w 


£- 

a 


1  1  1  1  *  1  •  •  •  1  i  1  1  1  1  1  #  1  1  1  1  1  1  1  i 

£,£,S!JS!,$®®®®®©©©®©®®®©©©©©©© 

coaoonnr-coco^Trnrr  —  cNoooocD^LD’r-inooLDLnn 
^®52r^®°"*nTTCN,!Tm’COJ0^00(£)OcNtrir5inco 

SnS£2iti!rl^,^S®r'nr'(71"‘CNCNnronn^ 

<X>'r-’r-'r-CiCVCr)nCr)Cr)’*TitTt'<T*rTtrtTT-'T'V'<T*trT*TTf 


CO 


■a 

b* 

•H 

tn 

«J 

w 

u 

<D 

+J 

0) 


u 

d 

ft 


CO 


«3 

0 

(ft 

I 

L 

a 


11 

£ 

+-• 

a 

® 

T3 


© 

U 

(ft 

C 

J 

1 

c. 

a 


CO 

CM 

CM 

CN 

1 

© 

) 

® 

© 

« 

© 

00 

CT) 

00 

ao 

CO 

ID 

00 

r*' 

00 

Tf 

in 

r- 

00 

r- 

O)  00 

c* 

CM 

in 

a> 

i  co  <j>  cm 


LT)incOtO(£>(jC>U3(0(X)(X)tP(O(f) 


oflSSS?!NP,CNtNCNNtNNtNNNWWCNNN«Ni 


©©©©©©©©©©©©a)©©®©® 


JSSSSfS1?®®®®®00  00  ®00®®®®®®®®®® 

®®(O(OCO(OiO(D(p(OU)(I>(O(PiP(j}(O(O(Oi£}(p(0^}^)(^ 


a 

0 


~  88888888888888888881 
2  8888888888888888888555555 


lOOOOO 

iinoiflOin 

>  O  —  CN  CN 


> 

o 


72 


ON-AXIS  INTENSITY 


ltfqr 


Fig.  27:  Normalized  intensity  on  the  beam  axis  as  function  of  detector 
FOV  for  the  same  parameters  as  in  Fig.  20  except  RD  =  0.2  cm. 
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o  d) 


RECEIVED  POWER 


FOV  (RAD) 

Fig.  28:  Received  power  (arbitrary  units)  on  the  beam  axis  as  a  function  of 
detector  FOV,  for  the  same  parameters  as  in  Fig.  27.  The  received 
power  is  multiplied  by  exp  (x) . 
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TABLE  13.  Received  Power  and  On- Axis  Intensity  as  Functions  of  the  Detector  FOV  for  the 
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ON-AX  13  INTENSITY 


Fig.  29:  Normalized  intensity  on  the  beam  axis  as  a  function  of  detector  FOV 
for  the  same  parameters  as  in  Fig.  20,  except  RD  =  2.0  cm. 
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RECEIVED  POWER 


Fig.  30:  Received  power  (arbitrary  units)  on  the  beam  axis  as  a  function 
of  detector  FOV,  for  the  same  parameters  as  in  Fig.  29.  The 
received  power  is  multiplied  by  exp  (t) . 
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4.  CONCLUDING  REMARKS 


In  this  report,  we  have  discussed  some  new  formulations,  computational 
methods,  and  results  for  the  complex  subject  of  multiple  scattering  effects 
on  laser  beams  traversing  dense  aerosols,  particularly  in  connection  with 
topics  such  as  higher  order  scattering  contributions,  exact  solutions  of  the 
radiative  transfer  equation  in  the  small  angle  approximation,  and  scattered 
radiance  as  a  function  of  detector  f ield-of-view.  Conclusions  based  on  work 
done  thus  far  and  recommendations  for  future  research  are  given  as  follows. 

1.  In  this  paper,  we  developed  the  formulations  for  solving  the 
radiative  transfer  equation  (RTE)  for  laser  beams  traversing  scattering  media 
in  the  small-angle  approximation  and  for  making  the  results  tractable  to 
numerical  computations  by  three  different  approaches: 

a.  Tam  and  Zardecki  exact  approach 

b.  Box  and  Deepak  exact  approach 

c.  Fante-Dolin  approximation  approach 

The  formulations  for  the  quantities  of  interest,  namely,  the  beam  intensity, 
irradiance  and  the  detected  power  are  given.  Numerical  computations  for  the 
scattered  irradiance  and  detected  power  were  performed  using  the  first  two 
approaches.  Tam  and  Zandecki  approach,  though  computationally  more  time 
consuming  than  the  Box  and  Deepak  approach,  gives  additional  insight  into 
the  contributions  made  by  higher  orders  of  scattering  to  the  total  multiple 
scattering;  the  latter  just  gives  the  total  MS  value.  To  date,  the  Tam 
and  Zardecki  approach  has  not  been  used  to  solve  the  case  of  non-Gaussian 
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phase  functions;  whereas  the  Box  and  Deepak  approach  can  handle  non-Gaussian 
phase  functions.  It  is  recommended  that  further  work  be  performed  to 
solve  for  the  radiance  and  irradiance  of  laser  beams  for  non-Gaussian 
phase  functions  by  these  two  approaches. 

2.  It  should  be  noted  that  the  small-angle  approximation  is  valid  for 
highly  forward-peaked  phase  functions.  It  is,  therefore,  recommended  that 
approaches  be  developed  to  deal  with  the  case  of  MS  in  laser  beams  traversing 
small  size  aerosol  particles,  with  broadly-peaked  phase  functions. 

3.  It  is  recommended  that  experiments  and  numerical  computations  be 
performed  to  determine  the  effects  of  the  detector  field  of  view  on  the 
extinction  and  backscattering  function  measurements. 

4.  It  is  recommended  that  measurements  and  theoretical  computations  be 
performed  to  investigate  the  effects  of  beam  diameter  on  optical  extinction 
measurements . 

5.  It  is  recommended  that  the  above  investigation  be  repeated  for  the 
case  of  backscattering,  which  is  of  great  importance  for  single-ended 
electro-optical  systems  such  as  Lidar,  DIAL  and  Laser  Doppler  Velocimeter  (LDV) . 
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ATTN:  STEEP-MM-I  S 

ATTN:  STEEP-MT-DS  (CPT  Decker) 

Ft.  Huachuca,  AZ  85613 

Commander 

Dugway  Proving  Ground 

ATTN:  STEDP-MT  (Dr.  L.  Solomon) 

Dugway,  UT  84022 

DEPARTMENT  OF  THE  NAVY 
Commander 

Naval  Research  Laboratory 

ATTN:  Code  5709  (Mr.  W.  E.  Howell) 

ATTN:  Code  6532  (Mr.  Curcio) 

ATTN:  Code  6532  (Dr.  Trusty) 

ATTN:  Code  6530-2  (Mr.  Gordon  Stamm) 

ATTN:  Code  8320  (Dr.  Lothar  Ruhnke) 

ATTN:  Code  8326  (Dr.  James  Fitzgerald) 

ATTN*  Code  43202  (Dr.  Hermann  Gerber) 
4555  Overlook  Avenue,  SW 
Washington,  DC  20375 

Chief,  Bureau  of  Medicine  &  Surgery 
Department  of  the  Navy 
ATTN:  MED  3C33 
Washington,  DC  20372 

Commander 

Naval  Air  Systems  Command 

ATTN:  Code  A  I R— 30 1 C  (Dr.  H.  Rosenwasser) 

Washington,  DC  20361 

Commander 

Naval  Sea  Systems  Command 

ATTN:  SEA-62Y 1 3  ( LCDR  Richard  Gilbert) 

ATTN:  SEA-62Y21  (A.  Kanterman) 

ATTN:  SEA-62Y21  (LCDR  W.  Major) 

Washington,  DC  20362 

Project  Manager 

Theatre  Nuclear  Warfare  Project  Office 
ATTN:  TN-09C 

Navy  Department 
Washington,  DC  20360 

Institute  for  Defense  Analysis 

400  Army-Navy  Drive 

ATTN :  L.  B i berma n 

ATTN:  R.  E.  Roberts 

Ar  1  i ngton ,  VA  22202 


Commander 

Naval  Surface  Weapons  Center 

Dahlgren  Laboratory 

ATTN:  DX-21 

ATTN:  Mr.  R.  L.  Hudson 

ATTN:  F- 56  (Mr.  Douglas  Marker) 

Dah I gren,  VA  22448 

Commander 

Naval  Intelligence  Support  Center 
ATTN:  Code  434  (H.  P.  St.Aubin) 

4301  Su  it  land  Road 
Suitland,  MD  20390 

Commander 

Naval  Explosive  Ordnance  Disposal 
Technology  Center 
ATTN:  AC-3 

Indian  Head,  MD  20640 

Of  f i cer- i n -Charge 
Marine  Corps  Detachment 
Naval  Explosive  Ordnance  Disposal 
Technology  Center 
Indian  Head,  MD  20640 

Commander 

Naval  Air  Development  Center 

ATTN:  Code  2012  (Dr.  Robert  Helmbold) 

Warminster.  PA  18974 

Commander 

Naval  Weapons  Center 

ATTN:  Code  382  (L.  A.  Mathews) 

ATTN:  Code  3882  (Dr.  C.  E.  Dinerman) 

ATTN:  Code  3918  (Dr.  Alex  Shlanta) 

China  Lake,  CA  93555 

Commanding  Officer 
Naval  Weapons  Support  Center 
Applied  Sciences  Department 
ATTN:  Code  50C,  Bldg  190 

ATTN:  Code  502  (Carl  Lohkamp) 

Crane,  IN  47522 

US  MARINE  CORPS 

Commanding  General 
Marine  Corps  Development  and 
Education  Command 
ATTN:  Fire  Power  Division,  D09 1 

Quantico,  VA  22134 
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DEPARTMENT  OF  THE  AIR  FORCE 


HQ  AFLC/LOWMM  1 

Wr ight-Patterson  AFB,  OH  45433 

HQ  AFSC/SDZ  1 

ATTN:  CPT  D.  Riediger 

Andrews  AFB,  MD  20334 

USAF  TAWC/THL  1 

Egl In  AFB,  FL  32542 

USAF  SC 

ATTN:  AD/YQ  (Dr.  A.  Vaslloff)  I 

ATTN:  AD/YQO  (MAJ  Owens)  1 

Eg  I  In  AFB,  FL  32542 

AFAMRL/TS 

ATTN:  COL  Johnson  1 

Wr I ght-Patterson  AFB,  OH  45433 


Comma  n  der 

Hanscom  Ai r  Force  Base 

ATTN:  AFG L-POA  (Dr.  Frederick  Volz)  1 

Bedford,  MA  01731 

Headq  ua rters 
Tactical  Air  Command 

ATTN:  DRP  1 

Langley  AFB,  VA  23665 

AFOSR/NE 

ATTN:  MAJ  H.  Wlnsor  1 

Bo  I  I Ing  AFB,  DC  20332 

Dr.  Charles  Arpke  1 

OSV  Field  Office 

P.0.  Box  1925 

Egl In  AFB,  FL  32542 

OUTSIDE  AGENCIES 

Battel le,  Columbus  Laboratories 
ATTN:  TACTEC  1 

505  King  Avenue 
Co  I umbu s ,  OH  4320 1 


Dr.  W.  Michael  Farmer,  Assoc. Prof.,  Physics 
University  of  Tennessee  Space  Institute  1 
Tul lahoma,  TN  37388 

ADDITIONAL  ADDRESSEES 

Office  of  Missile  Electronic  Warfare 
ATTN:  DELEW-M-T-AC  (Ms  Arthur)  1 

White  Sands  Missile  Range,  NM  88002 

US  Army  Mobility  Equipment  Research  and 
Development  Center 

ATTN:  DROME-RT  (Mr.  0.  F.  Kezer)  1 

Fort  Be!  voir,  VA  22.060 


Di rector 

US  Night  Vision  and  E0  Laboratories 
ATTN:  DRSEL-NV-Vl  (Dr.  R.  G.  Buser)  1 

ATTN:  DRSEL-NV-Vl  (Mr.  R.  Bergemann)  1 

ATTN r  DELNV-V I  (Luanne  Obert)  1 

ATTN:  DELNV-L  (D.  N.  Spector)  1 

Fort  Belvolr,  VA  22060 


Comma  nder 

217th  Chemical  Detachment 

ATTN:  AFVL-CD  1 

Fort  Knox,  KY  40121 

Headquarters 

US  Army  Medical  Research  and 
Development  Command 

ATTN:  SGRD-RMS  1 

Fort  Detrlck,  MD  21701 

Commander 

US  Army  Environmental  Hygiene  Agency 
ATTN:  Librarian,  Bldg  2100  1 

Aberdeen  Proving  Ground,  MD  21010 

Commandant 

Academy  of  Health  Sciences,  US  Army 
ATTN:  HSHA-CDH  1 

ATTN:  HSHA-IPM  2 

Fort  Sam  Houston,  TX  78234 


Toxicology  Information  Center,  JH  652 
National  Research  Council  1 

2101  Constitution  Ave. ,  NW 
Washington,  DC  20418 
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